Numerical and Analytical Approach to the Quantum Dynamics of 
Two Coupled Spins in Bosonic Baths 

Alessandro Sergj*| 
School of Physics, University of KwaZulu- Natal, 
Pietermaritzburg, Private Bag X01 Scottsville, 
3209 Pietermaritzburg, South Africa 

Ilya SinayskiJ3 

Quantum Research Group, School of Physics, 
University of KwaZulu- Natal, Durban, 4001, South Africa 

Francesco Petruccionaj; 
Quantum Research Group, School of Physics and National Institute for Theoretical Physics, 
University of KwaZulu- Natal, Durban, 4001, South Africa 

Abstract 

The quantum dynamics of a spin chain interacting with multiple bosonic baths is described in 
a mixed Wigner-Heisenberg representation. The formalism is illustrated by simulating the time 
evolution of the reduced density matrix of two coupled spins, where each spin is also coupled to 
its own bath of harmonic oscillators. In order to prove the validity of the approach, an analytical 
solution in the Born-Markov approximation is found. The agreement between the two methods is 
shown. 
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I. INTRODUCTION 



For the sake of studying quantum information transport in solid state devices, the quan- 
tum dynamics of spin chains coupled to bosonic baths has attracted much attention in the 
recent scientific literature [l|, [2I, y|, Q, If], Q, Is]. Here, we show how a mixed Wigner- 
Heisenberg representation of quantum mechanics is particularly well-suited to the numerical 
simulation of such systems. This is illustrated by studying the time evolution of the re- 
duced density matrix of a minimal chain, composed of two spins, each coupled to a bath of 
harmonic oscillators. The temperature of each bath can be defined independently, so that 
nonequilibrium situations can be addressed with no further theoretical or computational 
efforts. The dynamics of the total systems, spins plus harmonic oscillators, is unitary and 
numerically exact. No Markovian or rotating waves approximations need to be invoked. 
Reduced operators are obtained simply by integrating the coordinates of the oscillators in 
Wigner phase space. Our numerical solution is compared with an analytical solution of the 
Markovian master equation of the two spins and good agreement is found. It is very easy 
to extend the algorithm to study longer chains and multiple bosonic baths. 

It is worth remarking that the mixed Wigner-Heisenberg representation that we adopt in 
this paper has been originally proposed for introducing a quantum-classical representation 
of systems immersed in gravitational fields and in plasma physics js]. In particular it has 
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12 1, and it 



been developed [lOj and applied to a variety of models in chemical physics 
has already been noted [13( that such a representation is exact in the case of (bosonic) bath 
of harmonic oscillators. 

This paper is organized as follows. Section [III illustrates the Wigner-Heisenberg rep- 
resentation of quantum mechanics. Section II III provides the details of the model we have 
studied. The numerical algorithms for the computer simulation is illustrated in Sec. IIVI The 
Born-Markov approximation for the master equation and details of the analytical solution 
are given in Sec. |V] Results of both our numerical and analytical studies are displayed in 
Sec. EH Finally, our conclusions are reported in Sec. IVII1 
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II. WIGNER-HEISENBERG REPRESENTATION OF QUANTUM MECHANICS 



Let us consider a system defined by the total Hamiltonian operator 



H = Hq + Hb + H; 



l SB 



where the subscripts S, B, and SB stand for subsystem, bath, and coupling, 



The Heisenberg equation of motion of the density matrix can be written as 14] 



(1) 



respectively 



d A i 
dt P = ~h 



H p 



B 



H 

P 



(2) 



where B is the antisymmetric constant matrix 



B 



1 
-1 



(3) 



Assuming that the bath Hamiltonian depends on a pair of canonically conjugated opera- 
tors X = (R, P), and the coupling has the form Hsb — Hsb{R), we can introduce a partial 
Wigner transform for the density matrix 



Pw(X) 



1 f^ { a-l m + i) 



(4) 



and for the generic bath-dependent operator x{R, P) 

X W (X) = J dze ip - z/h (R- 



z . . z, 

- 2 \x\R + - 2 ) 



(5) 



where X = (R, P) are canonically conjugated classical variables in phase space. Taking the 
partial Wigner transform of Eq. (j2J) 



Of 



pw{X,t) 



Hw{X) p 



■T>- 



Hw{X) 
pw(X,t) 



(6) 



where 



_ e f d tBxj d 







(7) 



In Equation ((7|) we have used the symbol d j = d / dXj to denote an operator of derivation 
(with respect to the phase space point coordinates) which acts on whatever stands on its 



3 



right. Analogously, d / acts on whatever stands on its left. Moreover, the summation over 
repeated indices must be performed in Eq. (JTj) and in the following. The mixed Wigner- 
Heisenberg form of the Hamiltonian operator, Hw, is 



Hw{X) = H s + H W:B (X) + H W:S b(R) ■ 



(8) 



Equation ([6]) provides a mixed Wigner-Heisenberg representation of quantum mechanics, 
where operators also depend on phase space (c-number) coordinates, which is completely 
equivalent to the usual Heisenberg representation. However, the difficulties associated to 
the solution of Eq. are formidable. Yet, for quadratic bath Hamiltonians 



N ( P 2 1 \ 



(9) 



where (Rj, Pi), I = 1, . . . , N, are the coordinates and momenta, respectively, of a system of 
N independent harmonic oscillators with frequencies u>j, and for interaction Hamiltonians 
of the type 

H w ,sb = V B (R) H' s , (10) 

where Vb{R) is at most a quadratic function of R and H' s acts only in the Hilbert space of 
the subsystem, Eq. can be rewritten using the antisymmetric operator matrix 



lin 



-1 



1 1 + fdjBjjdj 

fiiBvtj 



(11) 



Actually, it can be shown that for the class of Hamiltonians specified by Eqs. fl9J) and ( flO 



lin 



(12) 



holds exactly. For more general bath Hamiltonians and couplings, such a substitution 
amounts to performing a quantum-classical approximation [h]]. What matters here is that 
for the class of systems we are interested in the Eq. f|T2l is exact and provides via Eq. ([2]) a 
Wigner-Heisenberg formulation of quantum mechanics which can be numerically simulated 
employing algorithms previously developed within a chemical-physical context 15]. 
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III. MODEL SYSTEM 



The system we are interested in this paper is defined by the following subsystem Hamil- 
tonian 

Hs = -j**P&<? - ~ j z ^af (13) 

representing a chain of two quantum spins coupled to each other. The constants ji, with 

(k ) 

i = x,y, z, dictate the strength of the coupling between the spins. The operators a\ with 
% = x,y,z are the Pauli matrix operators for spin k s = 1,2. The bath Hamiltonian is 



H 



W,B 



2 N p2 
r I,k a 



EE 



(14) 



fc s =i/=i " 

The above Hamiltonian represents two independent harmonic oscillator baths with coordi- 
nates and momenta (Ri,k s , Pi,k s ), (where I — 1, N labels the oscillators and k s — 1, 2 labels 
the bath). The harmonic oscillator frequencies uj are taken to be bath- independent since 
we want to adopt two baths with identical spectral density. However, the baths can have 
different initial conditions (and eventually different temperature). The coupling is given by 



N 



r {ks) 



(15) 



Hw,sb = - E E c iRi,k s <? z 
k s =l 1=1 

showing that each spin is coupled to its own oscillator bath. 

The density matrix of the two-spin chains obeys the exact Wigner-Heisenberg equation 



0_ 

Of 



Pw 



h 



Hw Pw 



lin ' 



Hw 

Pw 



(16) 

The reduced density 



where H w = H s + H w ,b + H w ,sb is given by the sum of Eqs. 
matrix of the spin subsystem is given at all times by 

. 2 N 

Ps(t)= / II T[dX IM J w (X,t) . (17) 

For the calculation presented in this paper, we assume an initially uncorrelated density 
matrix, which, once partially Wigner transformed, takes the form 



where 



Pw(to) = Ps(t )p w ,B(X,t ) 



A ^ tanh(/^/2) 

fc s =l 7=1 



(18) 



7T 



x exp 



o tanh(/3 fcs o; J /2) 

-l *lw,B 



(19) 



and where H WB is defined in Eq. ( TT4")) and (3j s = (k B T kg ) 1 is the inverse temperature of 
each oscillator bath is the Boltzmann constant). 



IV. NUMERICAL ALGORITHM 

In cases in which the coupling Hamiltonian Hw,sb can be treated as a small perturbation 
(weak coupling), it is useful to represent the abstract Eq. ffl6|) in the adiabatic basis. Such 
a basis is defined by the eigenvalue equation 

(H s + H w ,SB)\a;R) = E a {R)\a; R) . (20) 

Hence, Eq. (1161) can be recast in propagator form 

Pw\^t)=Y.(e-^) a P ^\X)^ (21) 

0(3' ' M 

where 

i£aa' ,/3f3' = a , 5 aa ' 6/3/3' + Zia',/3/3' • (22) 

The operator iC° aa , is defined as 

iC° aa , = iuj aa > + iL aa > , (23) 
where u aa i = (E a (R) — E a >(R))/h and 

iL aa > = P-^ + l(F% + Fti).-^. (24) 



Fyy = —(a;R\dHw/dR\a]R) is the Hellmann-Feynman force 16]. The transition operator 
T a a',pi3' is purely off-diagonal and defined by 



P - daP [ 1 + 2 P-d a , -dP) W 



, 1 (E a , - E p ,)d* alp , d 
' P "W[ 1 + 2 P-d^, dp) S <* 



(25) 



where d a p = (a;R\ d /dR\(3\R) is the coupling vector between the adiabatic states \a;R) 
defined in Eq. ( |2"U]) . Above and in the following, the quantities are defined adopting scaled 
coordinates, according to the definition of the Hamiltonians in Eqs. f IT3lfT5l) . The operator 
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FIG. 1: Time evolution trace of the reduced density matrix element vs time (Pi = /?2 = 0.005). 
Initial density matrix p s (0) = |1,0)(1,0|. The error bars display the numerical error. 

Taa',/3/3' realises the quantum transitions of the subsystem due to the coupling to the bath. 
Assuming weak coupling, and for the sake of comparison to a Markovian master equation, the 
action of the transition operator will be disregarded: This amounts to perform an adiabatic 
approximation of the dynamical evolution of the spin subsystem. 

In the adiabatic approximation the evolution of the density matrix becomes simply 

p$(X,t) = e-^^(I) 

= e ~' fo d ™*«< ' e -*«w p <$ \X) . (26) 

Equation fl2"6"l) shows that using the adiabatic approximation, in the adiabatic basis, the 
evolution of the density matrix in the Wigner-Heisenberg representation can be calculated 
by propagating classical-like trajectories, under the action of the Liouville operator ( |24l ). 
and considering a phase factor integrated along the trajectory. The initial X coordinates, 
representing the quantum state of the bath in phase space, can be sampled from the initial 
density matrix ffT^j) . 

Although in the adiabatic approximation the dynamics is easily calculated in the adiabatic 
basis, for quantum information problems it is more convenient to consider the reduced 
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density matrix 

p%(t) = J dX-£U m (R)p$ '{X^CU-X^R) (27) 

aa' 

in the natural basis 1 1) = |1,1), \2) = |1,0), |3) = |0, 1), |4) = 1 , 0) . The matrix U appearing 
in Eq. (j27p is, of course, the rotation matrix from the adiabatic to the natural basis which 
can be constructed, as well known, by using the adiabatic eigenvectors as columns. 

Everything seems quite straightforward so far. However, the definition of U is somewhat 
arbitrary, since the columns can be evenly permuted, and the adiabatic eigenvectors in the 
Wigner-Heisenberg representation of quantum mechanics depend on the configuration point 



R. In addition, the LAPACK [17j] numerical routines, which we have used to calculate the 
eigenvectors, return a matrix U with the columns ordered corresponding to the increasing 
value of the eigenvalues. It turns out that this configuration-dependent permutation of the 
columns of U introduces fictitious dynamics, as can be verified by propagating the density 
matrix p = |1)(1|, defined in terms of the natural state "spin-up spin-up" of the spin chain, 
which should be left invariant under the action of the Hamiltonian H$ + Hw,sb, defined in 
Eqs. (H3D and (TISD. 

In order to solve this problem, it is sufficient to note that one would like to have a rotation 
matrix U as close as it could be to the matrix S formed by ordering the Cartesian basis 
vectors e? (in the present case j = 1, . . . , 4), with e 1 = [l ], e 2 = [ 1 0] and so 
on. Upon writing u Q for the adiabatic eigenvectors, a = 1, . . . , 4, one can define a metric 

g a,j = ( u a _ e j) . ( u a _ ^ _ ^ 

The definition of the metric in Eq. ( |28i) allows us to solve the ordering problem in a unique 
way. As a matter of fact, for each j, labelling the columns of the desired rotation matrix, we 
can look for the a which minimizes the metric g a ^\ This leads to the possibility of ordering 
the columns of U in such a way that this matrix is as close as it can be to £, and it effectively 
solves the numerical problem with the fictitious dynamics arising from the permutations of 
the adiabatic eigenvectors along the phase space trajectory. 

V. MASTER EQUATION FOR THE COUPLED SPINS 

A system with total Hamiltonian ([T]), obeying the Liouville (Heisenberg) equation of 
motion Q, can be studied in the weak coupling limit by performing the Born-Markov 
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FIG. 2: Time evolution of the reduced density matrix element vs time (f3± = 1, fa = 0.3). 
Initial density matrix rho s (0) = l^oX^ol with |^o) = (| 1, 1) — \l,Q))/sqrt2. The continuous line 
is the analytical solution. The filled circles display the results of the numerical calculation. 



approximation [18j. In such a case, the equation for the reduced density matrix becomes 



d8tl B [Hj&(t), M>(t - 8),%>{t) ® /3 B (0)]], 



(29) 



where the index J denotes the interaction picture with respect to the free Hamiltonians of 
the system and bath. The operator ps denotes the reduced density matrix of the system S 
and f>B is the density matrix of the reservoir B. 

After performing the rotating wave approximation over the rapidly oscillating term in 
the master equation one gets: 

j t Ps(t) = -i[H s ,Ps(t)] + 

1 



(30) 



To obtain Eq. (I30p one assumes that the system-environment interaction has the form 
Hsb = Hi V% ® fi\ the operators Vi = if} and /j = j\ act on the system and the bath 
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FIG. 3: Time evolution of the reduced density matrix element vs time = fh = 0.005). 
Initial density matrix p s (0) = |1,0)(1,0|. The continuous line is the analytical solution. The filled 
circles display the results of the numerical calculation. They are joined by a dashed line to help 
the eye. 

degrees of freedom, respectively. In Eq. (1301) a Lamb-type renormalization Hamiltonian was 
neglected and decay rates 7 a ,^(iv) are given by the Fourier image of the bath correlation 
functions: 

7«*(") = l +C ° dse^{ft{s)f^)). (31) 

J — oo 

The transition operators V a (u) originates from the decomposition of the operator V a in 
the basis of the eigenoperators of the system Hamiltonian Hg. If one denotes the eigenvalues 
of the Hamiltonian Hg by e and the corresponding projection operator as 11(e) then: 

v a (u)= ]T n(e)y a n(e'). (32) 

e'—e=u> 

To obtain the master equation for the open system we rewrite the Hamiltonian of the 
whole system in the following way: 

H = H s + Hbi + Hb2 + Hsbi + H SB 2, (33) 

where ^5 is defined in Eq. ( TT3"]) and here we further assume that j x = j y = j, so that the 
constants j > and j z > denote the strenght of XY and ZZ interaction, respectively. 
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As already stated in the previous section, in this article scaled units are chosen, so that 
he — U — 1. We rewrite the Hamiltonians of the reservoirs k s — 1, 2 as 

H B k s = Y. u n,k s b{ M b ntks . (34) 

n 

The interaction between the spin subsystem and the bosonic baths is described by 



H S B ks = -^^g^(b n>k3 + blJ. 



(*.)/ 



To derive an equation of the form (I30p for the Hamiltonian 
eigenvalues and eigenvectors of the Hamiltonian H$ ffT3"]) : 



namely 



Hs = Ai|Aj)(Aj 



|Ai) = |1, 1), Ai = -j z , 

|A 2 ) = |0,0), X 2 = -j z , 

|A 3 ) = -^(|1,0) + |0,1)), A 3 = -2j+j 2 

|A 4 ) = -^(-|1,0> + |0,1», \ 4 = 2j+j 2 



(35) 

one needs to find the 
(36) 



(37) 
(38) 
(39) 

(40) 



In this basis, the transition operators take the form: 

|Ai)(A 1 |-|A 2 )(A 2 | 



yO) _ t/(2) 

Vq — Vq 



with u>o = 0. The operators Vq cause decoherence. The operators V 

V(i) = — |A 3 )(A 4 |, 
V(2) = |A 3 )(A 4 |, 



(41) 

(42) 
(43) 



describe the dissipation between the levels A 3 and A 4 with the transition frequency uj = 4j. 
Finally, the master equation takes the form: 



^ = -i[H s ,p} + j2(£ Di (p)+£ci(p)) 



i=l 



where 



t 1 h>t 



+ 



(v; t) pv (t) - ^ [v (i) v; t) , p 



(44) 

(45) 
(46) 



11 



and 



Caip) = (7 W (0+)+7 W (0-)) x 



1 



(47) 
(48) 



In the basis of eigenvectors of the Hamiltonian H$ the system of the corresponding 
differential equations can be solved. The exact solution of Eq. ( j44l) is 

i 



where we have introduced the elements of the density matrix: 

/n(t) = /u(0)> 



Mt) = / 22 (o), 



/»(*) 



1 -e 



+ 



-m\ /ss(O) 

/44(0) 



+ 



/ 44 (t) = (i - e" *) n_ 



n 



/as(0) 



-o/^ /44(0) 



+ 



/i 2 (t) = /i 2 (0) exp {-Ag c t + it{\ 2 - Ai)), 



(49) 



(50) 
(51) 

(52) 



(53) 



(54) 



/isW = /ia(0) exp - ft + t + it(\ 3 - AO), 



/u(t) = /u(0) exp - + zf(A 4 - Ai)), 



/ 23 (0 = / 23 (0) exp (-</ c t - n + t + it{\ 3 - A 2 )), 



fm(t) = / 24 (0) exp (-<? c t - Q_t + zi(A 4 - A 2 )), 



fu{t) = / 34 (0) exp (-nt + it{\ A - A 3 )). 



(55) 
(56) 
(57) 
(58) 
(59) 
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In the above expressions we have defined the constants 



± 2 



1 -(j( 1 \±u) + ^(±u)), (60) 



n = n + + n^, (6i) 

S c =^E(7 W (0 + )+7 W (0-)) • (62) 



2f =1 

The above solution will be used in the following as a reference for the numerical simulation. 

VI. CALCULATIONS AND RESULTS 

We have performed various numerical calculations varying the temperatures of the oscil- 
lator baths and compared to the analytical solution given in Sec. |VJ The coupling constants 
in the Hamiltonian (jT3"j) have been taken as j x = j y = j = 1 and j z = 1/2. The two baths, 
with iV = 200 harmonic oscillators each, have been assigned an Ohmic spectral density. To 
this end we employed the form of the coupling constants c/ and frequencies ojj introduced 
in Ref. ljj: 

ci = (^ujoujj) 1 ^ (63) 
u; 7 = - In (1 - Iu ) (64) 

where u = (1 —exp(—u max )/N, with £ = 0.007 and a; max = 3. In order to compare with the 
analytical solutions of the weak-coupling master equation of Sec. |V] we have performed an 
adiabatic propagation in the mixed Wigner-Heisenberg representation of quantum mechanics 
and sampled 50000 initial conditions to calculate the reduced density matrix, ps of the two 
coupled spins. Figure [1] shows the numerical precision of our numerical scheme displaying 
the constancy of the trace of ps versus time in the case of /?i = /3 2 = 0.005. 

In general, we have found a very good agreement between the results provided by both 
the numerical and the analytical approach for all the various temperatures investigated. 
Here, we discuss explicitly two calculations. 

Calculation (i) has been performed with the baths in a nonequilibrium configuration, at 
the two different temperatures f3\ = 0.3 and @2 = 1. The initial reduced density matrix ps{0) 
has been taken equal to ps(0) = |\l/o)(^o|, with |\l/ ) = ^(|1,1) — |1,0)). Figure [5] shows 
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the comparison between the numerical and the analytical dynamics of the matrix element 
pf. In this case, the analytical solution is 



Of course, in such a low-temperature case the Markovian approximation is expected to 
provide very good results and this is numerically confirmed. 

Calculation (ii) has been performed with fli — 02 — 0.005 and an initial p s equal to 
p s = |1,0)(1,0|. Figure [3] shows the comparison between the numerical and the analytical 
dynamics of the matrix element p| 2 . The theoretical solution is in this case 



At higher bath temperature, the Markovian approximation (used in the analytical solution) 
can describe the numerical results in a good but qualitative way. The difference in the 
oscillation frequencies of the analytical and the numerical solutions arises from neglecting the 
Lamb-type renormalization of the Hamiltonian Hs in the derivation of the master equation 
in the Born-Markov approximation. The discrepancy in the long time decay of the analytical 
and numerical results arises from the fact that 02, to) in the analytical expression of 

p 22 (t) should contain some memory effects on the time-interval on which the evolution is 
considered. 

VII. CONCLUSIONS 

Upon adopting a mixed Wigner-Heisenberg representation, we have shown how the quan- 
tum dynamics of two coupled spins interacting with multiple bosonic baths can be numer- 
ically simulated. An analytical solution in the Born-Markov approximation has also been 
found and we have shown agreement between these two approaches. 

Both the analytical and the numerical method can be generalized in order to study 
additional coupled spins, in order to build longer spin chains immersed in independent 
bosonic baths. Equilibrium and nonequilibrium situation can be addressed on an equal 
basis. 

The numerical algorithm is suited to include nonadiabatic correction in the unitary evo- 
lution of the density matrix of the total systems. As such, it can also be used to assess novel 
approaches to non-Markovian dynamics of open quantum systems. 




(65) 
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